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Abstract 

Tomato late blight caused by the oomycete pathogen Phytophthora infestans (Mont.) de Bary is a major threat to tomato 
production in cool and wet environments. Intensified outbreaks of late blight have been observed globally from the 1980s, 
and are associated with migration of new and more aggressive populations of P. infestans in the field. The objective of this 
study was to reassess late blight resistance in the wild tomato accession L3708 {Solanum pimpinellifolium L.) against 
pathogens of different aggressiveness. An F23 genetic mapping population was developed using L3708 as the paternal 
parent. Two isolates of P. infestans, Pi39A and Pi733, were used for inoculation. Pi733 is a highly aggressive genotype that 
defeats three known late blight resistance genes, Ph-1, Ph-2, and Ph-5t in tomato. In contrast, Pi39A is a less aggressive 
genotype that defeats only Ph-1. Restriction site Associated DNA Sequencing (RAD-Seq) technology was used to massively 
sequence 90 bp nucleotides adjacent to both sides of PstI restriction enzyme cutting sites in the genome for all individuals 
in the genetic mapping population. The RAD-seq data were used to construct a genetic linkage map containing 440 single 
nucleotide polymorphism markers. Quantitative trait locus (QTL) analysis identified a new disease-resistant QTL specific to 
Pi733 on chromosome 2. The Ph-3 gene located on chromosome 9 could be detected whichever isolates were used. This 
study demonstrated the feasibility and efficiency of RAD-Seq technology for conducting a QTL mapping experiment using 
an F2:3 mapping population, which allowed the identification of a new late blight resistant QTL in tomato. 
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Introduction 

Late blight, caused by the oomycete pathogen Phytophthora 
infestans (Mont.) de Bary, is a devastating disease affecting tomato 
{Solanum lycopersicum L.) and potato production especially in cool 
and wet environments. Intensified epidemic outbreaks of tlie 
disease have occurred throughout the world since the 1 980s. This 
is associated with the migration of new and more aggressive 
pathogen populations [1,2]. The global predominant genotype of 
P. infestans in both tomato and potato field before the 1980s was 
designated as a single lineage, US-1 [3]. Newer immigrated 
genotypes of P. infestans isolates are usually highly aggressive and 
resistant to metalaxyl fungicides; they quickly displace the original 
US-1 genotype [3,4]. Dramatic population shift of P. infestans was 
occurred in Taiwan from 1998, and the highly aggressive isolate of 
the US-1 1 clonal lineage displaced the original US-1 clonal lineage 
[5]. The two isolates, Pi39A and Pi733 used in this study represent 



the population shift of P. infestans in Taiwan. Pi39A belongs to the 
US-1 clonal lineage [5], whereas Pi733 belongs to the US- 11 
clonal lineage. They were a part of collections at AVRDC to 
survey P. infestans populations in Taiwan from 1997 to 2008. 

Genetic factors associated with resistance to late blight in 
tomato have been characterized in several wild tomato species 
[6,7,8,9,10,11]. Among these resistance genes, Ph-3 has been 
widely used in tomato breeding programs as it confers resistance to 
P. infestans isolates in many regions [12,13]. Ph-3 was originally 
identified from the wild tomato Solanum pimpinellifolium accession 
L3708 and maps to the distal end of chromosome 9, close to the 
DNA marker TG591 [8,14]. Two studies using advanced tomato 
lines derived from L3708, implied that in addition to Ph-3 at least 
one other gene contributes to late blight resistance in L3708. The 
first study found that advanced lines containing the resistant Ph-3 
allele of L3708 were overcome in the field, but wild-type L3708 
plants were not [15]. The second study demonstrated that one 
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group of advanced lines conferred stronger resistance against 
highly aggressive P. infestans isolates than a different group of 
advanced lines, even though all lines had the same homozygous 
Ph-3 genotype [12]. Therefore, it is important to determine 
whether there is a new genetic factor associated with resistance to 
late blight in L3708. 

Quantitative trait locus (QTL) mapping has long been the 
standard method to identify resistance genes of wild tomato 
accessions against late blight [7,8,9,10,11]. Despite the develop- 
ment of genetic mapping populations and the measurement of 
phenotypes, typical QTL mapping requires great effort to identify 
new polymorphic markers and their genotypes, for all individuals 
in a mapping population, when new genetic crosses are made 
[8,16,17,18,19,20,21,22,23,24,25]. Restriction site Associated 
DNA sequencing (RAD-seq) [26] is a new sequencing-based 
genotyping method. It is able to overcome difficulties in identifying 
polymorphic markers, especially for crosses between accessions 
with low genetic polymorphism. The technical basis behind RAD- 
seq technology is to sequence a few million short DNA sequence 
reads anchoring at specific restriction enzyme sites in the genome. 
Alignment of sequencing reads between two parental genotypes 
allows detection of polymorphic single nucleotide polymorphism 
(SNP) sites as DNA markers, and identification of marker 
genotypes for all individuals in a mapping population, with or 
without the reference genome sequences [27]. Furthermore, RAD- 
seq reduces costs through multiplexing of bar-coded individuals 
[26]. RAD-seq technology has been successfully applied in 
numerous studies in crops requiring construction of genetic maps 
and QTL analysis [28,29,30,31,32,33,34]. 

This study reexamines the genetic components of resistance to 
late blight resistance in wild tomato S. pimpinelltfolium L3708 by 
QTL mapping using RAD-seq technology. Our results demon- 
strated the feasibility and efficiency of RAD-Seq technology in 
conducting a QTL mapping experiment using an F2:3 mapping 
population in crops with reference genomic sequences, and 
identified a new late blight resistant QTL on chromosome 2 of 
L3708. 

Materials and Methods 

Plant materials 

An F2:3 genetic mapping population was developed from the 
cross of S. pimpindlifolium accession L3708 (resistant parent) with 
the elite cultivar S. lycopersicum T3224 (susceptible parent). Four 
hundred and sixty nine individual F2 plants were grown in the field 
of the experimental farm of Known-You Seed Co LTD in Tainan, 
Taiwan for collection of F3 seeds. DNA ^\as c-xtracted from 
individual Fj plants for RAD-seq. One hundred and twenty F2:3 
Unes were randomly chosen to assess disease reaction by 
inoculation of two different P. infestans isolates. Six accessions, 
TS33 (S. lycopersicum L6161), TS19 (S. lycopersicum L6160), L3708 
{S. pimpinellifolium L3708), LA1()33 {S. hahrochaites VI60()17), 
WVa700 {S. lycopersicum var. cerasiforme L6193) and CLN2037B 
(S. lycopersicum AVTO9808) were used as differential hosts to 
identify physiological race of the two P. infestans isolates [5]. The 
tomato cultivar M82 was used as a susceptible control to assess 
disease severity for individual lines of the genetic mapping 
population [35]. 

Seeds were sown individually into 45 mm diameter plastic pots 
containing 1:1 volume mixture of peat-Ute (King Root Plant 
Medium #3, Taiwan) and peat moss. Liquid fertilizer 
(N:P:K= 30:10:10, HYPONEX #5, USA) was applied every 
other week. Seedlings were raised in a room at the Phytotron of 
the National Taiwan University at 25/20°C (day/night) for four 



weeks. They were then moved to a growth chamber for 
inoculation and symptom development. 

Assessment of disease reaction to late blight 

P. infestans isolates Pi39A and Pi733 were used for the 
inoculation studies. Pi39A belongs to the US-1 clonal lineage 
and was re-isolated from Pi39-inoculated tomato plants at 
AVRDC [5]. Pi39 was originally collected in Tainan, Taiwan in 
1997. In contrast, Pi733 belongs to the US-11 clonal lineage - a 
newly immigrated genotype, and was collected in Nantou, Taiwan 
in 2007. They were a part of collections at AVRDC to survey P. 
infestans populations in Taiwan from 1997 to 2008. The isolates 
were maintained on rye A agar plates at 18°C [36]. Inoculum was 
prepared from mycelial cultures grown on rye A agar plates at 16- 
18°C in the dark for 12-14 days. A sporangial suspension diluted 
to 5 X 10* sporangia per milliliter was incubated at 12°C for 2 h to 
induce zoospore release. Plants were inoculated using a sprayer to 
atomize the zoospore/ sporangia suspension onto the foliage to the 
point of run-off, and incubated at 20°C with 100% relative 
humidity in the dark for 24 h. Plants were then maintained at 60- 
95% relative humidity, with a daily 14 h light period. 

Symptoms of late blight was visually scored based on a disease 
severity rating (DSR) of 0-6: where 0 indicates no symptoms; 1 
indicates <5% leaf area affected and small lesions; 2 indicates 6- 
15% leaf area affected and restricted lesions; 3 indicates 16-30% 
leaf area affected and/or water-soaked flecks on stems; 4 indicates 
31-60% leaf area affected and/or a few stem lesions; 5 indicates 
61-90% leaf area affected and expanding stem lesions; 6 indicates 
91-100% of leaf area affected, extensive stem damage, or a dead 
plant [8] . Plants were individually scored for DSR when the M82 
susceptible control displayed the most severe symptoms. 

To determine physiological race of each P. infestans isolate, 24 
seeds from each of the 6 differential hosts and the susceptible 
control M82 were sowed, inoculated and evaluated together. One 
inoculation was carried out for each isolate. The mean of DSR 
score of a differential host was used to determine disease reactions 
against the inoculated P. infestans isolate. Disease reaction is 
designated as "resistant" when the mean DSR is less than 3. In 
contrast, disease reaction is designated as "susceptible" when the 
mean DSR is equal or largc-r than 3. 

For the QTL mapping experiments, 6 plants of each of the 120 
F2:3 lines, and 18 plants of the susceptible control M82 were 
evaluated together at each inoculation; four inoculation trials were 
carried out for each isolate. The mean DSR score of an F2:3 line 
inoculated with the same isolate was the average of the mean DSR 
scores from 4 trials and was used for data analysis in the QTL 
analysis. However, if less than 4 plants of an individual F2:3 line 
were germinated in a trial, no mean DSR score was calculated. 
The overall mean DSR score of an individual F2:3 line was treated 
as a missing value on the occasion that the mean DSR score were 
calculated less than 3 times out of 4 repeat trials. 

Construction of the RAD library 

Total genomic DNA was isolated from young tomato leaves 
using a modified CTAB method [37], and further purified using 
the DNeasy Blood & Tissue Kit (QIAGEN, Venlo, Netheriand) 
following manufacturer's instructions. 

At/-digested RAD libraries were prepared following the 
protocol of Etter et al. [38]. Sixty F2 samples were multiplexed 
with 60 different PI barcodes in one RAD hbrary (Table SI). For 
each sample, 1 |Xg gDNA was digested with 20 units of Pstl-HF 
(New England BioLabs [NEB], Ipswish MA, USA) overnight in a 
50 \xL reaction volume. Samples were heat-inactivated for 20 min 
at 80°C. Digested DNAs were ligated to 2 JiL 100 nM PI bar- 
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coded adapters, a modified Solexa adapter (Table SI), along with 
1 HL lOx NEBufier4 (NEB, Ipswish MA, USA), 0.5 |xL 2000 unit 
HL"' T4 DNA ligase (NEB, Ipswish \L\, USA), and 0.6 nL 
100 niM rihoATP (Promega, Madison WI, USA) in a 60 |J,L 
reaction volume for 1 h. Samples were heat-inactivated for 20 min 
at 65°C. Samples were pooled (20 nL each, 60 samples) and a 
50 [lL aliquot was loaded into a 0.5 mL PGR tube (Axygen 
catalog # PCR-05-C, Corning, Tewksbury, MA, USA). DNAs 
were .sheared using a Bioruptor UCD-200 sonicator (Diagenode, 
Liege, Belgium) set to high, for 3 run.s of 7 min (30 s on/30 s off). 
The peaks of most DNA aliquots were approximate 300 bp. If the 
peak of sheared DNA was over 500 bp, then additional sonications 
were performed until the peak became less than 500 bp. Sheared 
DNA aliquots were pooled and concentrated using two MinElute 
columns (QIAGEN, Venlo, Netherland) and eluted with 40 nL EB 
buffer (10 niM Tris-Gl, pH8.5) in each column. Eluted DNAs 
were combined and size selected using Agencourt AMPure XP 
magnetic beads (Beckman Coulter, Brea CA, USA) with a volume 
DNAibeads ratio of 1 :0.65; this removed DNA fragments of less 
than 300 bp [39]. Recovered DNAs were suspended in 20 ^L EB 
buffer, and treated using a Quick Blunting Kit (NEB, Ipswish MA, 
USA) for end repair. 1 \lL Blunt Enzyme Mix, 2.5 JtL 10 x 
Blunting Bufier, and 2.5 \lL 1 mM dNTP mix were added to the 
20 |xL DNA solution. The mixture was incubated at 25°C for 
30 min. Agencourt AMPure XP magnetic beads (Beckman 
Coulter, Brea GA, USA) were used for reaction clean-ups with a 
volume DNAibeads ratio of 1 : 1 .8; this removed DNA fragments of 
less than 50 bp. The repaired dsDNAs were suspended in 20 \lL 
buffer EB and quantified using the Quant-iT dsDNA HS Assay 
Kit (Life Technologies, Carlsbad CA, USA). Adenine was added 
to the 3' ends of dsDNA fragments in a 50 )xL reaction volume 
containing 1 |xg dsDNAs, 5 |J,L 10 x NEBuffer2, 1 ixL 10 mM 
dATP, and 3 (xL of 5 unit |a,L~ ' Klenow Fragment (NEB, Ipswish 
MA, USA) mbced and incubated at 37°G for 30 min. DNAs were 
cleaned using 90 [iL (1.8 X volume) Agencourt AMPure XP 
magnetic beads (Beckman Coulter, Brea CA, USA), and 
suspended in 45 |xL EB buffer. Reactions for P2 adapter ligation 
were assembled by adding 1 |j.L 10 |J.M P2 adapter (Table SI) to 
the dsDNA solution along with 5 ^lL 10 x NEBufier2, 0.5 ^L 
100 mM riboATP (Promega, Madison WI, USA), and 0.5 jxE 
2000 unit |xL"' T4 DNA ligase (NEB, Ipswish MA, USA). The 
mixture was incubated at 20°C for 3 h. The P2 adapter-ligated 
dsDNAs were then purified using 35 |j,L (0.7 x volume) Agencourt 
AMPure XP magnetic beads (Beckman Coulter, Brea CA, USA), 
suspended in 20 |a,L EB buffer, and quantified using a Quant-iT 
dsDNA HS Assay Kit (Life Technologies, Carlsbad CA, USA). 
50 ng of this DNA product was PCR amplified using 4 |j.L 10 |J,M 
modified Solexa primer mix (Table SI) and 50 |a,L Phusion High- 
Fidelity PCR Master Mix (NEB, Ipswish MA, USA) in a 100 nL 
reaction volume. The PCR setting was 98°C for 30 s, followed by 
18 cycles of 98°C for 10 s, 66°C for 30 s, 72°C for 30 s, and a final 
extension reaction at 72°C for 5 min. The PCR-enriched product 
was purified with 70 |jL (0.7 x volume) Agencourt AMPure XP 
magnetic beads (Beckman Coulter, Brea CA, USA), and 
normalized to 10 ng ^iL One RAD library was sequenced in 
one lane of an lUumina Hiseq2000 flow cell (100 bp single-end 
reads) (lUumina Inc., San Diego, CA, USA). Next-generation 
sequencing (NGS) was provided by the Genome Research Center 
at the National Yang-Ming University, Taiwan. 41 1,738,303 reads 
were obtained (Table S2). Sequences are available at the Sequence 
Read Archive http://www.ncbi.nlm.nih.gov/Traces/sra/, at ac- 
cession SRAl 44571 (Table S2). 



Sequence analysis 

Stacks vl.08 (http://creskolab.uoregon.edu/stacks/) [27] and 
CLG Genomics Workbenc h v6.5.1 (CLC Bio, http://www.clcbio. 
com) were used for sequence analysis. Raw sequencing reads were 
processed using the "process_radtags" command in Stacks to filter 
out reads with quality scores less than Q,10, and to sort reads to 
individual samples based on barcode sequences. Sorted reads for 
each sample were aligned to the tomato genome sequence buUd 
SL2.40 (SOL Genomics Network, http://solgenomics.net/), using 
the read mapper tool in the CLC Genomics Workbench. Stringent 
parameters were used to j)rc'vc'nt high false positive rates for SNP 
calling. As most RAD read sequences had low possibility of 
overlap, the parameters for sequence alignment were set to allow 
no more than two mismatches for the 90 bp short reads (length 
fraction = 1.0, similarity fraction =0.97), and to discard aligning 
results if reads were mapped to more than two positions in the 
genome. The other parameter settings for read mapping were 
mismatch cost = 2, insertion cost = 3, and deletion cost = 3. The 
sequencing-read-alignment files of the two parents were used to 
define the SNP sites, and those of the 120 progenies used to 
determine genotypes at the defined SNP sites. This was performed 
using the "ref_map.pl" command of the Stacks software, set at 
default parameters, except for minimum read depth, which was set 
to 3. 

Development of additional SNP markers for genotyping 
by VeraCode technology 

In order to obtain independent genotypic data different from 
the RAD-seq data to delimit the chromosomal region where 
recombination events were suppressed, 1 14 SNP sites between the 
two parental lines were successfully developed as a customized 
genotyping panel using VeraCode technology (Illumina Inc., San 
Diego, CA, USA) (Table S3). At beginning, a total of 144 
candidate SNP sites were chosen from a SNP database that was 
independendy buUt from the RAD-seq data. The SNP database 
contains approximately 0.7 million SNPs distributed evenly in the 
tomato genome and was obtained by a comparison between two 
assembled genomic sequences of the parental Unes, S. pimpinelli- 
folium L3708 and S. lycopersicum T3224, each of which was 
subjected to whole genome shotgun sequencing with approxi- 
mately 8x coverage. To select the candidate SNP sites for 
Veracode SNP markers, two adjacent SNP sites were chosen with 
an approximately physical distance of 6 Mb on each of the 12 
tomato chromosomes. Genotypes of the 144 VeraCode SNP 
markers for individuals in the F2:3 mapping population were 
determined using the GoldenGate Genotyping Assay on the 
BeadXpress Reader System (Illumina Inc., San Diego, CA, USA) 
[40]. Fluorescence signal data were analyzixi using GenomeStudio 
genotyping module vl.O (Illumina Inc., San Diego, CA, USA). 
Sixteen markers failed to generate genotypes, and a further 14 
markers showed no polymorphism in the mapping populations. 
Reproducibility of genotypes of the VeraCode SNP markers can 
reach to 99.9%. It hence provide great accuracy to delimit the 
recombination suppression regions. A substantial number of RAD 
markers located in these chromosomal regions were able to be 
ignored, so it greatiy enhanced the efficiency to run the JoinMap 
software which is used to construct the genetic linkage map. 

Genetic map construction, and QTL analysis 

JointMap v4. 1 software was used to construct genetic maps 
[41]. Regression mapping was used as the mapping algorithm. A 
linkage group was defined when markers showed recombinant 
frequency smaller than 0.4, and independence LOD values larger 
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than 7.0. The Haldaiie mapping function was used to calculate 
genetic distance [42]. 

MapQTL v6 software was used for QTL analysis [43]. Interval 
mapping was first used to identify markers that significantly 
associated with phenotypic variations in the mapping population. 
Next, the MQM mapping module in the MapQTL software 
(similar to composite interval mapping) was used to confirm the 
final result of the QTL analysis by iterating significant cofactor 
(marker) selection. A regression algorithm was used to calculate the 
approximate LOD scores for interval mapping, MQM mapping, 
and the permutation test. The mapping step size was set to 
1.0 cM. The empirical threshold of LOD score to claim QTLs was 
set to 4.0; this was obtained by running tiie permutation test 
50,000 times. 

Results 

Designation of physiological races of Phytophthora 
infestans isolates 

Pi733 (US- 11 clonal lineage) and Pi39A (US-1 clonal lineage) 
isolates were used in this study to investigate the resistance profile 
of L3708 against different clonal lineages of P. infestans. Designa- 
tion of physiological race for these two isolates was based on the 
disease reaction of six differential tomato hosts: TS19 {ph+), TS33 
WVa700 {Ph-2), CLN2037B {Ph-3), L3708 {Ph-3,4i), and 
LA1033 {Ph-5t) [5]. Pi39A only infected TS19 and TS33, and was 
designated as race Tl. Pi733 infected TS19, TS33, WVa700, and 
LA1033, and designated as race Tl,2,5t (Table 1). 

Phenotypic evaluation of disease reaction for the 
mapping population 

To identify resistance QTLs to different P. infestans isolates, an 
F2:3 genetic mapping population was developed from a cross of 
resistant S. pimpinellifolium accession L3708 and susceptible 
cultivated tomato S. lycopersicum. One hundred and twenty Fj::? 
lines were randomly selected for phenotypic evaluation. Seedlings 
grown from the 120 Fj:^ lines were inoculated separately with 
either Pi39A or Pi733 isolates. Disease severity ratings for each 
inoculation were scored at the time that the susceptible tomato 
cultivar M82 control reached maximum severity. The DSR scores 
for resistant parent L3708 were 0.112 and 0.144, for Pi39A and 
Pi733, respectively. The DSR scores of the 120 ^2:'.', lines against 



Pi39A ranged from 0.109 to 4.984, with a mean of 1.537, and 
skewed towards zero (Figure lA; Table S4). In contrast, DSR 
scores of the 120 ^2-'i lines against Pi733 ranged from 0.346 to 
5.921, with a mean of 3.149, and showed an approximate equal 
frequency among the different DSR classes (Figure lA; Table S4). 
These results implied that the Pi733 isolate was more aggressive 
than the Pi39A isolate. The correlation coefficient between the 
DSR scores when inoculating Pi39A and Pi733 was 0.704 
(Figure IB), which implied that common genes conferred 
resistance against both isolates. 

Genotyping and construction of genetic map 

The total read number of the 120 bar-coded F2 samples 
obtained from two lanes of the lUumina Hiseq2000 platform was 
396,100,974. Approximate 92% of reads passed the default quality 
filter in the Stack software, and 76% of reads aligned to the tomato 
genome sequence build SL2.40. The aligned read numbers of 
individual samples in the F2 population ranged from 842,457 to 
4,446,219, with an average of 2,504,121 (Table S2). The RAD 
reads of the two parental samples were generated from an 
additional RAD library. There were 5,925,802 and 6,384,883 
aligned reads for the resistant parent S. pimpinellifolium accession 
L3708 and the susceptible parent S. lycopersicum, respectively. A 
total of 67,339 distinct 90 bp DNA sequences on the tomato 
genome sequence build SL2.40 were aligned with RAD reads 
from both parents, and 12,718 sites of 90 bp sequences were 
defined as RAD markers showing distinct homozygous haplotypes 
between the two parents. The genotypes of these RAD markers for 
each of the 120 Fj samples were determined, however, a 
substantial portion of genotypes for all RAD markers were missed. 
To remove potentially problematic RAD markers from use for 
constructing the genetic map, three artificial criteria were set: (I) 
markers with more than 1 0 missing data points were removed; (2) 
for the goodness-of-fit test of unbiased genotype segregating ratio, 
markers with a chi-square value larger than 30 were removed; (3) 
markers not mapping with reference sequences of the 12 
chromosomes were removed. After filtering, 4697 RAD markers 
remained. For ease of identification, the RAD markers were 
designated using a "OOgOOOOOOOO" format. The first two digits 
indicate the chromosome number, while the last eight digits 
indicate the physical position of the PstI cutting site on tomato 
genome reference sequence build SL2.40. 



Table 1. Race designation for Phytophthora infestans isolates Pi39A and Pi733. 
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"DSR: disease severity rating. 
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doi:l 0.1 371/journal.pone.009641 7.t001 
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Figure 1. Disease severity rating (DSR) for isolates Pi39A and 
Pi733. (A) Distribution of DSR in the mapping population. (B) 
Relationship of the two DSR score sets generated from independent 
inoculations of Pi39A and Pi733. 
doi:10.1371/journal.pone.0096417.g001 

The nature of RAD markers tends to produce low proportions 
of miscalled genotypes; these mainly result from insufficient read 
depths of a defined RAD marker [44]. To mitigate genotyping 
error, 24 samples were removed using the number of missing 
genotype as the selection criterion. Most removed samples also 
contained a lower number of RAD reads (Figure 2). Genotypes of 
the remaining 96 F2 samples were used to construct the genetic 
map. Additional 94 genotypes belonging to customized VeraCode 
SNP markers were added to the genotypic data. This allowed 
precise delimitation of chromosomal regions where very low 
recombination occurred for each of the 12 linkage groups. Among 
the 4697 RAD markers, 2047 contained no more than 2 missing 
genotypes out of the 96; these were selected for linkage analysis, 
along with the 94 VeraCode SNP markers, using JoinMap 
software. From the linkage analysis, the 2141 markers were 
separated into 12 linkage groups using a LOD threshold 
parameter of 7.0. Markers with physical locations identified on 
the same chromosome were grouped into the same linkage group, 
except for marker 06g33932831, which was grouped with markers 
on chromosome 9. The cause of this exception remains to be 
determined. 

To construct a concise accurate hnkage map, only a part of the 
RAD markers were selected. The selection criteria included: (1) 
allow the order of markers on the genetic map to be the same as 
the order on the physical map; and (2) the genetic distance 
between adjacent markers was in the range of 2 to 5 cM if a 
suitable marker was available. To achieve both criteria systemat- 
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Figure 2. Relationship between read number and missing data. 

Circles indicate each line of the F2 mapping population. The 96 filled 
circles represent samples used for genetic map construction, while the 
24 open circles represent samples removed because of a high number 
of missing data in the population. 
doi:10.1 371/journal.pone.009641 7.g002 

ically, RAD markers from each linkage group were selected by the 
following strategy. Markers were ordered by physical position on 
the color-coded genotype panel in the JointMap program. Possible 
genotyping errors were identified by looking for individual 
genotypes different from neighboring marker genotypes. Cross- 
over events could be identified when one genotype the same for 
more than three adjacent markers along the marker order 
switched to a different genotype identical for the following three 
markers in the same individual. The genotyping errors of each 
marker were counted, and crossover events along the marker 
order were recorded for the whole mapping population. In 
recombination blocks defined by four to seven crossover events, 
the marker with the lowest number of possible miscalling 
genotypes was selected. 

The final genetic map for the 12 chromosomes included 395 
RAD markers and 45 VeraCode SNP makers (Figure 3; Table S5 
and S6). The average genetic distance between two adjacent 
markers was 3.56 cM. None of the adjacent marker pairs had 
genetic distances larger than 13 cM. This indicated that the 
marker density of the whole genome was sufficient to capture 
major genetic effects causing phenotypic variations for QTL 
analysis. 

QTL analysis 

To understand whether the wild tomato accession L3708 
confers additional genes resistant to highly aggressive isolates of P. 
infestans, phenotypic data collected from inoculation of Pi39A and 
Pi733 were used separately for QTL analysis. Because the same 
genetic materials and genotypic data were shared for different 
QTL analyses, comparison between two QTL analyses could 
reveal distinct genetic factors contributing to new resistance 
against the highly aggressive isolates of P. infestans Pi733. 

A QTL at the distal end of chromosome 9 was detected when 
the Pi39A isolate was used (Figure 3; Table 2 and S7). The 
location of this QTL is the same as the Ph-3 gene [8] . The Ph-3 
QTL explained 41.5% of the phenotypic variation in the mapping 
population. The additive effect was —0.78, and the dominant 
effect was —0.26. This indicated that the L3708 allele increased 
resistance to the Pi39A isolate and led to a reduction in the DSR 
score when the L3708 allele replaced the allele from cultivated 
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Figure 3. Genetic map and results of QTL analyses for resistance to tomato late blight. The naming of the RAD markers was in the format 
"OOgOOOOOOOO". The first two digits indicate chromosome number, while the last eight digits indicate the physical position of the PstI cutting site on 
the tomato genome reference sequence build SL2.40. The naming of the VeraCode SNP markers was designated in the same way as the RAD markers, 
except for replacing "g" with "n", and the last eight digits indicating the physical position of the SNP site. VeraCode SNP markers names are 
highlighted green. The red bar indicates the location of the QTL resistant to isolate Pi733, and the hollow bar indicates the location of the QTL 
resistant to isolate Pi39A. 
doi:l 0.1 371/journal.pone.009641 7.g003 



tomato. The L3708 allele is partially dominant to the allele from 
the cultivated tomato. 

Two QTLs were detected when the disease reaction of the 
aggressive isolate Pi733 was examined (Figure 3, Table 2 and S7). 
The QTL on chromosome 9 mapped to the same chromosomal 
region as that to the Pi39A isolate, so this QTL was also 
designated as the Ph-3 gene. This Ph-3 QTL explained 63. 3"/!) of 
the phenotypic variation. The additive effect was —1.70 and the 
dominant effect was —0.18. This indicated that the L3708 allele of 
Ph-3 increased resistance to the isolate Pi733, and the effect of the 
L3708 allele of Ph-3 was almost additive when L3708 allele 
replaced the allele from cultivated tomato. The second detected 
QTL was close to the proximal end of the centromere on 
chromosomal 2. This QTL was not identified in the first round of 
interval mapping analysis, but detected after the marker closely 
linked to the Ph-3 gene was used as a cofactor in the MQM 
mapping analysis. This QTL has not been previously reported and 
was named qPh2.1. qPh2.1 explained only 1 1.9% of the phenotypic: 
variation. The additive effect was —0.63 and the dominant effect 
was —0.31. Therefore, the L3708 allele of qPh2.1 was partially 
dominant to the cultivated tomato allele of qPh2.1. The L3708 
allele of qPh2.1 also confers resistance to the Pi733 isolate. When 
compared with the effect oiPh-3 in the same mapping population, 
qPh2.1 can be classified as a minor QTL. 

Discussion 

The RAD linkage genetic map 

Two factors are major keys for SNP discovery using next 
generation sequencing technology. The first is to build correct 

alignment between reads and the reference sequences, the second 
is to generate a sufficient number of read sequences from the high- 
throughput sequencing platform [45] . The former can reduce the 
false positive rate, while the latter can reduce the false negative 
rate. For the RAD protocol, correct alignment is easier to achieve; 
this is because all RAD sequence reads have the same length (90bp 
after trimming barcode sequences and restriction site remains) and 
anchored at specific restriction sites [PstI), thus greatly reducing 
alignment ambiguity. Furthermore, correct SNP sites show normal 
genotypic segregation in the progenies. Therefore, it is easy to fix 
incorrect alignment problems using the RAD data. In contrast, an 
insufficient number of short read sequences resulting from the 
nature of RAD library constructions or lack of restriction sites in 
one parental line, is the major causes of missing data and 
genotyping errors [44]. 

To obtain a sufficient number of short read sequences for a 
sample in a cost-effective manner, four factors were considered 
before RAD fibraries were generated: (1) the number of PstI 
restriction sites in the tomato genome; (2) the read depth for a 
given RAD site; (3) the number of barcoded samples pooled in a 
RAD library; and (4) the product of the aforementioned three 
factors was set equal to the total number of single-end short reads 
generated from one lane of the Illumina Hiseq2000 platform. We 
performed a conservative estimation for these factors as follows: (1) 
There are 82,814 P.stI sites in tomato genome sequence build 
SL2.40, therefore 165,628 RAD sites were expected; (2) the 
average read depth for all RAD sites was set to nine, based on an 



estimation of 99yo probability that the read depth of any RAD site is 
no less than three; (3) 60 barcoded tomato samples were pook-d in a 
RAD library; (4) at least 89.4 miUion reads per lane were required to 
be generated from the Illumina Hiseq2000 platform, and each 
sample required at least 1.5 miUion reads in order to obtain correct 
marker genotypes. In this study, the Illumina service center was able 
to generate 198 million reads per lane. However, 12 samples had 
the read number less than 1.5 million. These samples also had at 
least 327 missing genotypes out of 4697 markers (Figure 2). To 
mitigate genotyping error problems [44], we set 327 missing 
genot)'pes as the threshold, and removed additional 12 samples 
whose number of missing genotypes greater than 327. Therefore, 96 
samples were remained to construct the genetic linkage map. 

The length of the genetic map in this study is longer than in 
previous studies [25,46]. Inffation of genetic maps usually results 
from genotyping errors [47]. The same phenomena was observed 
for construction of an Arahidopsis genetic map in a backcross 
population using genotyping technology similar to RAD sequenc- 
ing [30]. Because marker order can be affected by genotyping 
errors, we deployed a strategy to identify markers that resulted in 
the same marker order between the genetic and physical maps. 
Consequentially, the overall genotyping error rate for the 395 
RAD markers was reduced to 2. 18'M), a percentage that would not 
seriously affect the QTL analysis results. 

The genetic map in this study showed similar features to a 
typical tomato genetic map in which high recombination rate was 
found at the distal ends of all chromosomes, but recombination 
was suppressed in the pericentromeric regions [25,48]. AU 
VeraCode SNP markers in the pericentromeric regions showed 
identical genotypes for all individuals in the mapping population, 
therefore, markers on the region with identical genotype could be 
defined as the region that recombination was suppressed. 
Furthermore, the physical regions at which recombination was 
suppressed were similar to regions found in a previous study (Table 
S3) [25]. 

It is worth to note that our strategy to select RAD markers 
showing coUinear orders between the physical map and the genetic 
map depends on how well the physical map was built. The higher 
accuracy the physical map was built, the more the coUinear RAD 
markers could be obtained. The predicted tomato genomic size is 
900 megabases (Mb), of which 760Mb were assembled into 91 
scaffolds with most gaps restricted in the pericentromeric region 
where recombination were commonly suppressed [48]. While 
some RAD sites might locate at gaps outside the scaffolds, these 
RAD markers were removed early by filtration of potentially 
problematic markers. In addition, base accuracy of the tomato 
reference genome sequences is approximately one substitution 
error per 29.4 kb [48]. Given that 82,814 Prf/ sites were identified 
in the tomato reference genome, only 507 base substitution errors 
were expected as the false positive markers which were counted 
less than one tenth of overall genuine SNPs. Finally, a recent study 
which constructed a high density genetic map using the tomato 
high-throughput genotyping array with 7720 SNPs, demonstrated 
good collinearity between the genetic map and the physical map 
with few exceptional regions on chromosome 3, 10, and 12 [25]. 
The genetic map constructed in this study had medium marker 
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density and was obviously not affected by those exceptional 
chromosomal regions. In conclusion, the tomato physical map is 
well built, so our strategy to select coUinear RAD markers for 
genetic map construction is effective. 

QTLs associated with resistance to late blight 

The major QTL associated with resistance to late blight in 
L3708 was identified in the segment between 66,536,514 and 
67,494,653 bp on chromosome 9, with the most significant marker 
locating at 66,864,250 bp (Table 2). This QTL was assigned as the 
Ph-3 gene because TG591, the closely linked marker for the Ph-3 
gene, is located at genomic position 66.794 Mb (http:// 
solgenomics.net/) [8]. Furthermore, the Ph-3 gene has recendy 
been delimited within the chromosomal region from 66,714,091 to 
66,825,552 bp [49]. 

In this study, it seems that the Ph-3 gene is race non-specific, 
because the Ph-3 QTL was identified in both inoculations of Pi39A 
and Pi733. Nevertheless, from the result that examined the 
physiological race for Pi39A and Pi733, both isolates cannot 
infected the differential hosts L3708 and CLN2037B both of 
which possess the homozygous Ph-3 resistant allele (Table 1). 
Based on Flor's gene-for-gene theory [50], this result implied that 
both Pi39A and Pi733 isolates have the avirulent factor that is 
incompatible with the Ph-3 resistant allele. Identification of the Ph- 
3 QTL from inoculation of different isolates could result from the 
fact that the other parental line S. Ijcopersicum T3224 possesses the 
homozygous susceptible allele at the Ph-3 locus and hence the 
resistant allele and the susceptible allele at the Ph-3 locus were 
segregated in the ¥2-?, mapping population. Therefore, results in 
this study still support the conclusion in the previous study that the 
Ph-3 gene is a race-specific gene [8]. 

The main purpose of this study was to reexamine the genetic 
components of resistance to late blight in wild tomato S. 
pimpinellifolium L3708, and to ascertain whether another gene 
contributes to its late blight resistance. The results revealed that 
the newly identified qPh2.1 QTL confers resistance specifically 
against isolate Pi733. However, whether the qPh2.1 QTL is able to 
hold off the highly aggressive pathogen that can defeat the Ph-3 
gene [15], or if the qPh2.1 QTL merely enhances resistance against 
highly aggressive isolates, as shown in this and in Kim and 
Mutschler's work [12], remains to be fuUy determined. 

Supporting Information 

Table SI Information of barcoded adapters used to 
construct RAD libraries. 

pa.sx) 

Table S2 Total read numbers for each individual in the 
mapping population. 

(XLSX) 

Table S3 VeraCode marker information and genotypes. 
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Table S4 Disease reactions to P. infestans isolates 
Pi39A and Pi733. 
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Table S6 Genetic map information. 
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Table S7 QTL mapping results. 
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